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1.0  INTRODUCTION 


In  analyzing  a  physical  system  which  is  described  by  a 
system  of  differential  equations,  it  is  desirable  to  deter¬ 
mine  the  effects  that  variations  in  related  system  design 
parameters  have  on  the  system's  transient  response.  If  the 
governing  differential  equations  are  linear  with  constant 
coefficients,  then  an  investigation  of  the  effects  of  param¬ 
eter  variations  can  be  accomplished  by  examining  the  behavior 
of  the  roots  of  the  associated  characteristic  equation. 
Requiring  that  the  roots  have  negative  real  parts  is  both 
necessary  and  sufficient  for  stability  (Ref.  1) . 


For  example,  consider  a  physical  system  described  by 
the  following  equation: 

mx  +  2h k  +  k^x  =  0  (1) 

where  m,  h,  and  k  are  non-zero  design  parameters.  By  letting 

(2) 

(3) 

Equation  (1)  may  be  rewritten  as 


X1  “  x 
x2  -  *1  =  i 


X1  “  x2 


•  ••  2 

=  x  =  -  k  /m  x^  -  2h/m  x2 

which  can  be  expressed  as 


— 

— 

—  — 

“  «■ 

0 

1 

X1 

X1 

-2h 

X2 

x2 

m 

m 

(4) 

(5) 


(6) 


This  can  be  expressed  in  matrix  notation  as 

■> 

-►  • 

AX  =  X 


(7) 
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By  assuming  a  solution  of  the  form 


X  =  eXtY 

Eq.  (7)  becomes 


(8) 


(A-AI)Y  =  0  (9) 

where  I  is  the  identity  matrix.  Equation  (9)  represents  a 
system  of  linear  homogeneous  equations  and  has  a  non-trivial 
solution  if  the  determinant  of  (A-AI)  is  zero.  The  deter¬ 
minant  of  (A- XI)  is  referred  to  as  the  characteristic  equa¬ 
tion  and  is  given  by 

det(A-Xl)  =  X2  +  2h/m  X  +  k2/m  =  0  (10) 

If  the  roots  of  Eq.  (10)  have  negative  real  parts,  then  the 
system  described  by  Eq.  (1)  is  said  to  be  stable.  Hence  the 
solution  to  the  homogeneous  differential  equation  can  be 
determined  by  finding  the  eigenvalues  of  the  matrix  A.  If 
the  eigenvalues  have  negative  real  parts,  then  the  solution 
is  said  to  be  stable. 


The  coefficients  of  the  characteristic  equation  are 
functions  of  the  design  parameters  (m,  h,  and  k) .  The 
effects  of  variations  in  these  parameters  can  be  investi¬ 
gated  by  examining  the  behavior  of  the  roots  of  the  charac¬ 
teristic  equation.  This  process  is  generally  referred  to  as 
a  parametric  study. 


The  method  presented  herein  was  developed  to  investigate 
the  effects  that  variations  in  aerodynamic  parameters  have  on 
the  transient  response  characteristics  of  an  aircraft.  The 
method,  which  is  based  on  the  movement  of  the  roots  of  the 
characteristic  equation,  is  presented  in  generalized  form 
and  is  illustrated  with  simple  examples.  An  example  of  the 
technique  applied  to  investigating  aircraft  lateral  stability 
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is  also  presented.  For  variations  in  one  design  parameter, 
the  effects  on  stability  are  described  using  root  locus. 

For  simultaneous  parameter  variations,  the  effects  on  stabil¬ 
ity  are  represented  using  stability  region  maps. 

2.0  FORMULATIONS 


Consider  the  following  characteristic  equation  of 
degree,  n 

n 

G  =  l  F.<X)X1  =  0 
i=0  1 


where 


X  =  (x 


1 '  2  ‘ 


'V 


(11) 

(12) 


are  defined  as  design  parameters,  and  F^(X)  are  real  valued 
functions  which  have  continuous  first  partial  derivatives 
with  respect  to  X.  Note  that  X  is  a  function  of  X, 


Uk(X)}  =  £(X)  (13) 

and  is  chosen  such  that  G  is  identically  zero.  By  taking 
the  total  differential  of  G,  the  following  equation  is 
obtained: 


dG 


<!f>  axk  + 

K 


(|S_) 


where 


and 


+ 


-) 


,3G 

9  A  +m 
*k 


dxi  =  0 


3G  r 
TxT  =  l 


n  3F. 
^(X1) 
i=0  3xj 


1  <  j  <  l 


3G 

XX  - 


n  •  i  + 
l  iX1  1  F. (X) 
i=l 


(14) 

(15) 

(16) 
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As  indicated,  the  partial  derivatives  in  Eq.  (14)  are  to  be 
evaluated  at  a  point  (P“)  where 


<x?'x2' 


(17) 


For  a  characteristic  equation  of  degree  (n) ,  the  expression 
for  dG  represents  a  system  of  n  homogeneous  linear  equations. 
By  rewriting  Eq.  (14),  the  following  system  of  equations  is 
obtained: 


dX.  =  - 


dA  +m  l  dxa  pa  1 

*1  el 


dA~  =  - 


<lxy  <  <%>  +">+<i:>  dV 

Po  1  SL  £m 

e2  p2 


dX  =  - 
n 


(If K{<|2-)  dx,  +..  .+  (|S_)  dx.} 
n  n 


(18) 


which  can  be  expressed  in  matrix  form  as 


d t  =  AdX 


(19) 


where  dX  is  a  column  vector  of  n  components,  dX  is  a  column 
vector  of  l  components,  and  A  is  an  n  by  L  matrix  whose  ele¬ 
ments  are  complex.  In  general, 


Ai,j  -  - 

3  pj  pi 


(20) 


and  is  complex  whenever  X^  is  complex. 
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The  system  of  linear  equations  described  by  Eq.  (18) 
gives  a  linear  approximation  for  changes  in  each  root.  It 
is  valid  unless 


,36, 

W+k 

which  appears  in  the  denominator  of  the  Aj^j  element, 
vanishes.  It  may  be  shown  that 


l 

Wjk 


x 

vanishes  only  when  A ^  occurs  as  a  multiple  root  of  G. 


An  equivalent  way  of  obtaining  the  linear  approximation 
is  to  expand  the  function  (G)  in  a  Taylor  series  about  the 
point  (P  )  (Ref.  2) .  The  first  variation  of  G  gives  a  linear 
approximation  to  the  changes  in  the  roots  as  a  function  of 
changes  in  the  design  parameters.  The  resulting  linear 
system  is  again  given  by  Eq.  (19) . 

,  ,  ■>  -> 

By  approximating  dA  by  A A  and  dX  by  AX,  Eq.  (19)  becomes 

At  =  AAX  (21) 

By  specifying  AX,  the  corresponding  changes  in  AA  may  be 
determined.  Equation  (21)  represents  a  system  of  linear 
equations  which  relate  changes  in  each  root  as  a  function  of 
changes  in  the  design  parameters. 


The  characteristic  equation  (G)  is  an  algebraic  equation 
with  real  coefficients.  As  such,  there  are  definite  relation 
ships  between  the  coefficients  and  the  roots.  For  a  charac¬ 
teristic  equation  of  the  form  given  by  Eq.  (11) ,  the  product 
of  the  roots  is  given  by 

(-l)n  (A1A2  •••  An)  =  Fo(X)/Fn(X)  (22) 
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and  the  sum  of  the  roots  is  given  by 

-,*1  +  x2  +"-+  V  *  (23) 

Similarly,  other  relationships  exist  between  the  roots  and 
remaining  coefficients  of  G  (Ref.  3).  Either  Eq.  (22)  or 
Eq.  (23)  can  be  used  to  determine  the  error  in  the  roots 
introduced  by  using  the  approximation.  This  eliminates  the 
necessity  of  solving  the  characteristic  equation  to  estimate 
the  error  and  provides  a  means  whereby  the  approximating 
technique  given  by  Eq.  (21)  can  be  optimized.  The  error 
estimate  and  the  optimization  are  discussed  later  in  conjunc¬ 
tion  with  the  root  locus  and  boundary  mapping  applications. 

The  basis  of  the  root  movement  technique  is  the  system 
of  linear  approximating  equations  given  by  Eq.  (21) .  -  The 
effects  that  variations  in  the  design  parameters  have  on 
each  root  of  the  characteristic  equation  are  given  as  a 
weighted  sum  of  the  incremental  changes  in  the  design  param¬ 
eters,  where  the  weights  are  given  by  Eq.  (20).  The  rela¬ 
tive  weights  assigned  each  incremental  design  parameter  may 
be  used  to  rank  the  design  parameters  according  to  relative 
importance  for  each  root  at  a  particular  point.  For  systems 
containing  a  large  number  of  design  parameters,  the  ability 
to  select  only  the  most  significant  parameters  for  investi¬ 
gation  is  extremely  useful. 

3.0  MODIFIED  ROOT  LOCUS 

As  an  application  of  the  root  movement  technique,  con¬ 
sider  the  problem  of  determining  the  effect  of  variations 
in  one  design  parameter  on  a  system’s  transient  response 
characteristics.  The  effects  of  perturbations  in  one  de¬ 
sign  parameter  may  be  completely  described  by  examining  the 
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movement  of  the  roots  of  the  characteristic  equation. 

As  the  design  parameter  takes  on  values  ranging  from  plus 
infinity  to  minus  infinity,  the  roots  of  the  characteristic 
equation  move  in  the  complex  plane  along  trajectories  such 
that  the  characteristic  equation  remains  indentically  zero. 

A  complex  plane  plot  of  the  trajectories  of  the  roots  is 
referred  to  as  a  root  locus  plot.  By  examining  the  root 
locus,  the  effects  that  changes  in  a  particular  design  param¬ 
eter  have  on  both  the  frequency  and  damping  of  each  root  may 
be  described. 

The  process  of  .a  one-dimensional  parameter  investigation 
is  described  as  a  modified  root  locus  technique  to  indicate 
that  the  root  movement  approximation  described  in  Section  2.0 
is  used  to  generate  the  loci. 

To  initiate  the  process,  nominal  values  for  all  design 

parameters  are  selected,  and  then  the  characteristic  equation 

is  solved  to  determine  the  corresponding  roots.  That  is, 

-*-0  .  -*■  -y 

choose  X  and  find  A(X).  The  initial  starting  point  is 
referred  to  as  the  base  point,  and  the  corresponding  roots 
of  the  characteristic  equation  are  referred  to  as  the  poles. 
The  loci  emanate  from  the  poles  and  terminate  at  either 
infinity  or  at  points  in  the  complex  plane  referred  to  as 
zeros.  The  zero  points  are  defined  to  be  points  in  the  com¬ 
plex  plane  where  the  roots  of  the  characteristic  equation 
cluster  about  as  the  parameter  being  investigated  takes  on 
values  near  infinity.  For  a  characteristic  equation  of  the 
form  (G)  given  by  Eq.  (11) ,  the  cluster  points  or  zeros  are 
determined  as  the  roots  of  the  following  equation: 

n 

H  =  I 

i=0 


3F,  (X) 

5 1 

D 


X1  = 


(24) 
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where  Xj  is  the  design  parameter,  and  the  partial  derivatives 
are  to  be  evaluated  at  the  base  point. 


Starting  at  the  base  point,  the  incremental  changes  in 
each  root  are  calculated  using  Eq.  (21) ,  where  the  vector 
(AX)  contains  one  non-zero  entry  corresponding  to  the  design 
parameter  being  investigated.  For  example,  if  the  x^  design 
parameter  is  to  be  investigated,  then  AX  has  zero  entries 

everywhere  except  in  the  jth  position.  The  jth  position  of 

->■  ,  t 
AX  contains  an  entry  corresponding  to  the  incremental  change 

to  be  given  to  x j .  The  incremental  change  is  referred  to  as 

the  step  size.  As  noted  previously,  A  is  an  n  by  l  matrix 

with  complex  elements  of  the  form  given  by  Eq.  (20) ,  where 

the  partial  derivatives  are  to  be  evaluated  at  the  kth  point 

(Pk)  where 


(25) 


Notice  that  for  an  investigation  of  the  Xj  parameter, 

Pk  =  (x^x^,  •••  ,xk,  •••  ,x°,Xi)  (26) 

where  the  superscript  (0)  is  used  to  indicate  the  base  point 
values  of  the  other  l- 1  parameters.  By  using  the  incremental 
changes  found  above  for  the  roots,  the  roots  at  the  k  +  1 
point  are  defined  as 

tk+1  =  tk  +  Atk+1  (27) 

In  addition,  the  vector  (X)  at  the  k  +  1  point  is  defined  as 

Xk+1  =  Xk  +  A^k+1  (28) 
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Having  found  X  by  using  the  approximation,  the  error  in 
the  roots  may  be  estimated  as  follows.  The  relationships 
between  the  coefficients  of  G  and  the  roots  were  previously 
mentioned  in  Section  2.0.  If  G  is  an  nth  degree  polynominal 
of  the  form  given  by  Eq.  (11) ,  then  an  estimate  of  the  abso¬ 
lute  error  in  the  product  of  the  roots  at  the  kth  point  is 
defined  as  follows: 


error  = 


(-1) 


a  W 


vv 


(29) 


Similarly,  the  error  in  the  sum  of  the  roots  may  be  expressed 
as 


error 


A. +X-+* • *+A 

Vi  (xk> 

12  n 

F_  (X,. ) 

(30) 


If  the  error  at  the  kth  point  is  small,  the  above  process  is 
repeated  to  find  the  k+1  point.  However,  if  the  error  at 
the  kth  point  is  considered  excessive,  then  G  is  re-solved 
to  obtain  a  new  starting  point.  The  above  process  is  then 
repeated. 

It  is  reasonable  to  expect  that  a  variable  step  size  can 
be  used  in  determining  the  movement  of  the  roots.  With  the 
availability  of  an  error  estimate  at  each  step,  optimization 
of  the  technique  can  be  accomplished  as  follows.  If  the 
error  at  the  kth  point  is  larger  than  a  predetermined  toler¬ 
ance,  then  G  is  re-solved  to  establish  a  new  starting  point. 
If,  in  addition,  the  error  determined  at  the  k+1  point  is  also 
excessive,  then  the  step  size  is  halved  and  G  is  re-solved  to 
define  a  new  starting  point.  If,  on  the  other  hand,  the  error 
at  both  the  kth  and  k+1  points  is  less  than  a  predetermined 
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tolerance,  then  the  step  size  is  doubled  prior  to  finding  the 
k+2  point.  The  commitment  not  to  alter  step  size  unless  the 
same  error  condition  is  encountered  on  two  successive  steps 
is  intended  to  prevent  racing;  that  is,  constant  doubling  and 
halving. 

The  tolerance  used  in  determining  whether  or  not  the 
error  is  excessive  is  determined  empirically.  If  the  toler¬ 
ance  used  is  liberal,  then  the  trajectories  of  the  roots  will 
consist  of  linear  segments  which  are  piecewise  continuous. 
Decreasing  the  permitted  error  tends  to  smooth  the  approxi¬ 
mated  trajectories.  While  the  smoother  trajectories  are  a 
more  accurate  representation  of  the  loci,  the  increase  in 
computer  time  may  not  be  justifiable. 

As  stated  earlier,  the  loci  terminate  at  either  the 
zeros  or  at  infinity.  This  provides  a  convenient  way  to 
terminate  the  process.  It  is  only  necessary  to  require  that 
each  zero  point  contain  a  root  of  G  within  some  acceptable 
distance  of  it.  For  example,  the  process  can  terminate  at 
the  kth  point  only  if  each  zero  has  a  root  of  G  near  it. 

As  mentioned  in  Section  2.0,  the  system  of  approximating 
equations  is  invalid  at  multiple  root  points  because  the 
denominator  term  of  Eq.  (20)  vanishes.  In  the  complex  plane, 
the  multiple  root  points  are  seen  to  be  those  points  where 
the  roots  go  from  complex  to  real  or  vice  versa.  At  such 
points,  the  corresponding  tangent  vectors  have  infinite  slope, 
and  hence,  the  characteristic  equation  will  generally  require 
re-solving  at  least  once. 
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Example  I 


To  illustrate  some  of  the  ideas  discussed  concerning 
the  root  locus  technique,  consider  the  following  character¬ 
istic  equation: 

G  =  X2  +  F1(X)\  +  Fq(X)  =  0  (31) 

where  the  coefficient  functions  are 


and 


F1(X)=  xx  -  x2 
FQ(X)  =  0.25(x1  +  x2) 


(32) 

(33) 


where  x^  and  x2  represent  design  parameters.  Assume  that 
the  base  point  (P^)  is 


f° 


=  (xj,xfj) 


=  (3,2) 


At  P  ,  the  roots  of  G  are 


(34) 


X.  =  -  0.5  +  j 
and  1 

X2  =  -  0.5  -  j 


and  are  defined  to  be  the  poles  (P^  and  P2) .  The  zero 
points  corresponding  to  parameter  x^  are  determined  as 
the  roots  of  the  following  equation: 


(35) 

(36) 


1  3F.  (X)  . 

H  =  l  (-*± - )  X1  =  0 

i=0  3  $0 


=  ( 


3FQ  (X) 


"53? 


)  +  ( 


3FX(X) 


1  *0 


"53? 


1  ;o 


=  0.25  +  X 


(37) 
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Hence,  there  is  one  zero  point  at  X  =  -0.25.  The  loci  emanate 
from  the  poles  (P-^  and  P2)  and  travel  to  the  zero  and  to 
infinity.  The  trajectories  for  X^  and  X2  may  be  determined 
using  the  linear  approximation  given  by  Eq.  (21) .  For  this 
example,  Eq.  (21)  becomes 


OG/3xiy4G/8X)_)_0(3G/3x2  /8G/3X)^0  A^ 
P1  '  P1 


where 


(3G/3xiyW3X)^0(3G/3x2yW3X)_J,0  Ax2 


(xJ,x2,X1) 

(39) 

n- 

X 

M  O 

X 

ro  0 

to 

(40) 

aG/ax^^ 

■  X  +  0.25 

(41) 

3G/3x2 

=  -X  +  0.25 

(42) 

SG/SXj^  = 

2  \1  +  xx  -  x2 

(43) 

3G/3X2  = 

2  X2  +  X1  “  X2 

(44) 

At  the  base  point  (P  ) ,  Eq.  (38)  becomes 

f  AX.H  |"(0 . 5+0 . 125  j  )  (-0 . 5-0 . 375  j)l  [ax-j^ 


2  (0. 5-0. 125j)  (-0.5+0. 375j)  0 


Id 
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Since  Ax2  is  zero,  Eq.  (45)  simplifies  to 

AA^  =  (-  0.5  -  0.125j)Ax^  (46) 

A^2  =  (”  0.5  +  0.125j)Ax^  (47) 

In  the  complex  plane,  Eqs.  (46)  and  (47)  represent  vectors 
tangent  to  the  eigenvalue  trajectories  at  the  poles  P-j^  and 
?2*  letting  Ax^  equal  a  step  size  of  0.1,  an  approxima¬ 

tion  to  the  roots  at  P1  can  be  obtained  as  follows: 


t1  =  A0  +  At0 


t1 


(-  0.5  +  j) 

(-0.05  -  0. 012 j) 

+ 

(-  0.5  -  j) 

(-  0.05  +  0 . 012 j ) 

t1 


(-  0.55  +  0 . 99 j ) 
(-  0.55  -  0 . 99  j  ) 


(48) 


The  error  estimate  at  P1  may  be  computed  using  Eq.  (29)  as 

5.25 


error  = 


X1  A1 
A1  A2 


=  0.01 


(49) 


For  a  tolerance  of  0.1,  the  root  locus  shown  in  Fig.  1  was 
generated.  Table  1  summarizes  the  computations. 


The  root  locus  drawing  shown  in  Fig.  1  illustrates  the 
effects  of  increasing  parameter  (x^)  from  its  base  value  of 
3.0.  By  using  the  approximation  with  step  sizes  of  Ax^ =  0.1 
and  0.2,  it  was  necessary  only  to  re-compute  the  roots  once. 
This  occurred  in  the  vicinity  of  the  multiple  root  at  A =-1.28. 
In  sketching  the  approximated  trajectories  for  A;l  and  A2,  the 
circle  symbol  is  used  to  represent  the  approximated  A  value. 

The  vector  drawn  from  each  symbol  is  the  tangent  vector  com¬ 
puted  at  the  A  value  corresponding  to  that  symbol. 


17 


AEDC-TR-75-46 


By  using  a  negative  step  size,  the  root  locus  showing 
the  effects  of  decreasing  parameter  Xi  from  its  base  value 
of  3.0  may  be  generated.  Figure  2  shows  both  the  approxi¬ 
mated  and  actual  trajectories  for  this  case.  In  a  similar 
manner,  the  root  locus  showing  the  effects  of  increasing  and 
decreasing  parmater  X2  may  be  generated.  Table  2  summarizes 
the  computations. 


Figure  1. 


Root  locus  showing  the  effects 
of  increasing  design  parameter  x^ 
from  the  base  value. 


Figure  2.  Root  locus  showing  the  effects 

of  decreasing  design  parameter  xi 
from  the  base  value. 
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Table  1.  Summary  of  Computations  used  in  Obtaining  Figure  1 


X1 

*2 

Tangent 
Vector  (X^ 

Tangent 
Vector  (Xjl 

Approximate 

Roots 

Actual 

Roots 

3.0 

2.0 

-0.5  -0.125j 

-0.5  +0. 125 j 

- 

-0.5  -  0  j 

3.1 

2.0 

-0.5  -0.15 

3 

-0.5  +0.15  j 

-0.55  -0. 99j 

-0.55  ±0.953 

3.2 

2.0 

-0.5  -0.18 

3 

-0.5  +0.18  j 

-0.6  ±0.97j 

-0.6  -0.95j 

3.4 

2.0 

-0.5  -0.24 

3 

-0.5  +0.24 

3 

-0.7  -0. 93j 

-0.7  ±0.9  j 

3.7 

2.0 

-0.5  -0.35 

3 

-0.5  +0.35 

3 

-0.85  -0. 84 j 

-0.85  ±0. 85 j 

4.1 

2.0 

-0.5  -0.6 

3 

-0.5  +0.6 

3 

-1.05  -0.663 

-1.05  ±0.63j 

4.5 

2.0 

-0.5  -1.4 

3 

-0.5  +1.4 

3 

-1.25  -0. 36 j 

-1.25  ±0.25j 

4.7a 

2.0 

0.87+0 

3 

-1.87+0 

3 

-1.35  ±0.083 

-0.95,-1.75 

4.8 

2.0 

0.56+0 

3 

-1.56+0 

3 

-0.86,-1.94 

-0.9  ,-1.9 

4.9 

2.0 

0.5  +0 

3 

-1.5  +0 

3 

-0.8  -2.1 

-0.85,-2.05 

5.0 

2.0 

0.33+0 

3 

-1.3  +0 

3 

-0.75,-2.25 

-0.8  ,-2.2 

5.2 

2.0 

- 

- 

-0.68,-2.51 

-0.72,-2.48 

aRe-started  approximation  at  ■  4.7. 


Table  2.  Summary  of  Computations  used  in  Obtaining  Figure  2 


*1 

x2 

Tangent 
Vector  (Xj^) 

Tangent 
Vector (X2) 

Approximate 

Roots 

Actual 

Roots 

3.0 

2.0 

-0.5 

+0.125 j 

-0.5 

+0 . 125 j 

- 

-0.50 

±  0  3 

2.9 

2.0 

-0.5 

-0.1 

3 

-0.5 

+0.1 

3 

-0.45 

±1.01j 

-0.45 

±1.013 

2.8 

2.0 

-0.5 

-0.075j 

-0.5 

+0.075J 

-0.40 

-1. 02 j 

-0.40 

±1.023 

2.7 

2.0 

-0.5 

-0.05 

3 

-0.5 

+0.05 

3 

-0.35 

±1.033 

-0.  35 

±1.033 

2.5 

2.0 

-0.5 

-0 

3 

-0.5 

+0 

3 

-0.25 

±1.04j 

-0.25 

±1.033 

2.3 

2.0 

-0.5 

+0.05 

3 

-0.5 

-0.05 

3 

-0.15 

±1.  04  j 

-0.15 

±1.033 

2.1 

2.0 

-0.5 

+0.10 

3 

-0.5 

-0.10 

3 

-0.05 

±1. 03 j 

-0.05 

±1.013 

1.7 

2.0 

-0.5 

+0.2 

3 

-0.5 

-0.2 

3 

+0.15 

±0. 98 j 

+0.15 

±0.953 

1.5 

2.0 

-0.5 

+0. 253j 

-0.5 

-0. 253j 

+0.25 

±0. 94j 

•  +0.25 

±0.903 

1.3 

2.0 

-0.5 

+0.34 

3 

-0.5 

-0.34 

3 

+0.35 

-0. 89 j 

+0.35 

±0.843 

1.1 

2.0 

-0.5 

+0.43 

3 

-0.5 

-0.43 

3 

+0.45 

±0.823 

+0.45 

±0.763 

0.9a 

2.0 

-0.5 

+0.62 

3 

-0.5 

-0.62 

3 

+0.55 

±0.733 

+0.55 

±0.653 

0.5 

2.0 

-0.5 

+  1.4 

3 

-0.5 

-1.4 

3 

+0.75 

±0.363 

+0.75 

±0.253 

0.3a 

2.0 

+0.95+0 

3 

-1.95+0 

3 

+0.85 

±0.8  3 

+0.47, 

,  1.23 

0.2 

2.0 

+0.63+0 

3 

-1.63+0 

3 

+0.37, 

,+1.43 

+0.39, 

,+1.41 

0.1 

2.0 

- 

- 

+0.31, 

,+1.59 

+0.  33, 

,+1.57 

aRe-started  approximation. 
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4.0  STABILITY  REGION  MAPPING 


4.1  INTRODUCTION 

An  investigation  of  the  effects  of  variations  in  one 
design  parameter  is  illustrated  using  the  root  locus  tech¬ 
nique  described  in  Section  3.0.  The  effects  of  parameter 
variations  on  both  the  frequency  and  damping  may  be  obtained 
directly  from  the  root  locus  plot.  Since  it  is  conceivable 
that  the  loci  will  change  for  different  base  points,  general 
conclusions  concerning  these  effects  as  obtained  from  the 
root  locus  should  be  made  only  after  considerable  analysis. 

The  inability  to  generalize  is  a  major  disadvantage  of 
root  locus.  Stability  region  maps  provide  a  more  complete 
analysis  of  the  effects  of  parameter  variations.  For  the 
example  used  in  Section  3.0,  a  two-dimensional  stability  map 
completely  describes  the  effects  that  parameter  variations 
have  on  transient  response.  This  example  is  illustrated  at 
the  end  of  the  section. 

For  an  nth  degree  characteristic  equation  given  by 
Eq.  (11) ,  the  effects  of  simultaneous  variations  in  the  x^ 
and  Xj  parameters  on  stability  may  be  represented  by  a  sta¬ 
bility  region  map.  A  stability  region  in  the  x^,Xj  plane  is 
defined  to  be  the  set  points  (x^,Xj)  within  the  plane  which 
will  result  in  the  roots  of  the  characteristic  equation 
having  negative  real  parts.  The  region's  boundary  is 
defined  as  the  set  of  points  which  will  result  in  at  least 
one  root  of  G  having  a  zero  real  part  and  the  other  roots 
negative.  Stability  region  maps  serve  to  isolate  the  region 
of  stability  about  a  given  base  point. 

Stability  region  maps  may  be  refined  to  provide  addi¬ 
tional  information  concerning  the  roots.  Lines  of  constant 
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damping  may  be  constructed  within  the  region  using  the  root 
movement  technique.  In  addition,  information  concerning  the 
oscillatory  frequencies  is  available. 

A  stability  region  in  the  Xj ,x*  plane  encompassing  a 
■+0  ■L  J 

stable  base  point  (P  )  may  be  obtained  using  the  root  move¬ 
ment  approximation  given  by  Eq.  (21) .  The  method  used  in 
obtaining  the  stability  region  is  to  find  a  point  (P1)  on 
the  boundary  where  at  least  one  of  the  roots  of  G  is  zero 
and  the  other  roots  are  negative.  Once  a  point  on  the 
boundary  is  found,  the  approximation  is  used  to  find  other 
points  on  the  boundary  in  a  manner  such  that  the  boundary 
is  traversed. 

4.2  REGION  BOUNDARY 

The  first  point  on  the  boundary  (P1)  is  found  by  varying 
only  one  parameter.  If  x^^  represents  the  axis  of  abscissas, 
then  Axj  is  set  to  zero  and  Axj_  is  used  to  vary  x±  negatively 
from  its  base  point  value  until  the  boundary  or  a  predeter¬ 
mined  limit  is  found.  In  the  approximation  given  by  Eq.  (21)  , 
AX  contains  zeros  everywhere  except  in  the  ith  position. 

The  ith  position  of  AX  contains  the  negative  of  the  step  size. 
Beginning  at  the  base  point,  the  incremental  changes  in  each 
root  are  calculated  and  are  used  to  modify  X.  If  any  of  the 
roots  have  become  positive,  then  a  point  near  the  boundary 
has  been  found.  However,  if  the  new  roots  are  all  negative, 
then  the  approximation  is  used  to  determine  another  incre¬ 
mental  root  change.  The  above  process  is  repeated  until  a 
point  on  the  boundary  is  found. 

The  process  of  finding  the  first  point  on  the  boundary 

is  similar  to  the  modified  root  locus  technique  described  in 

Section  3.0.  The  only  difference  is  in  where  the  process  is 

.  . 

terminated.  As  in  the  root  locus  technique,  a  variable  step  size 
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is  used  where  the  step  size  is  adjusted  according  to  the 
criteria  established  in  Section  3.0.  The  first  point  on 
the  boundary  is  referred  to  as  P^,  where 


(50) 


where  the  superscript  (0)  refers  to  the  base  point  values. 


At  the  point  (P1) ,  assume  that  the  real  part  of  the  kth 
root  vanishes  and  that  the  other  n-1  roots  have  negative  real 


parts. 

The  linear 

-►1  . 

approximation  at  P  is  then 

AXX 

=—  (3G/3x^ y 

^3G/3X)  Ax.-(3G/3x.  /3G/3X)  Ax. 

-*■1  1  V  -*>1  3 

p  p 

*1  *1 

(51) 

ax2 

=-(3G/3xiy 

3G/3X)  Ax.  -  ( 3G/3x .  / 3G/3X)  Ax. 

*1  V  *i  3 

*2  *2 

(52) 

AX,  =-(3G/3x.  /3G/3X)  Ax.  -  (3G/3x .  /3G/3X)  Ax.  (53) 

K  x'  *1  1  3/  *1  3 


AX  =- ( 3G/3x.  / 3G/3X)  Ax.-(3G/3x.  /3G/3X)  Ax. 

n  1/  ji  1  V  £l  3 

n  n 

where 

£l  ,00  1  0  .  , 

Pk  =  (x1'x2'***'xi'’**'xji'X]c> 


(54) 


(55) 


-►1 

Since  Re(Xjc)  is  zero  at  P  ,  the  homogeneous  equation  for 

Re(AXk)  represents  a  vector  in  the  x^,Xj  plane  which  is  tan- 

->■1 

gent  to  the  boundary  at  P  .  By  moving  along  this  vector,  new 
points  on  the  boundary  may  be  found  where  Re(Xk)  remains  zero. 
The  objective  is  to  traverse  the  boundary  in  a  clockwise 
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manner  keeping  the  base  point  (P°)  to  the  right, 
this,  Eq.  (53)  is  rewritten  as 


To  achieve 


where 


al  = 


Re  C AXj^)  =  a^Ax^  +  a2AXj 
Im(AXk)  =  b^Ax^  +  b2AXj 


Re(Hr)Retlxl  +  Iir- <lf> 

Re(lf)  +  lm{^} 


(56) 

(57) 


(58) 


a~  = 


JL 


Re<lf>  + 


(59) 


and 


bi  ■ 


Be(|f)  +  lm(||) 


Im<lf7>Re<lf>  - 

b2 - 1 - — T- 


Re  (■|y^  (■|y) 


(60) 


(61) 


where  all  partial  derivatives  are  evaluated  at  P^.  From  the 

1  .  7  K 

point  (P  ) ,  a  second  point  on  the  boundary  (P  )  may  be  found 

by  moving  along  the  tangent  vector  given  above  by  Eq.  (56) . 

-H  ->-2  . 

From  P  ,  P  is  determined  by  letting  either  Axi  or  AXj  in 
Eq.  (56)  equal  the  step  size  and  then  solving  for  the  other 
Ax  required  to  satisfy  the  homogeneous  equation. 


The  slope  of  the  tangent  line  is  used  to  determine  which 
Ax  term  is  to  be  solved  for.  If  in  the  xi,Xj  plane,  x^ 
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represents  the  axis  of  abscissas,  then  for  tangent  lines  whose 
slope  magnitude  is  greater  than  one,  the  Axj  term  should  be 
adjusted  by  the  step  size  and  the  Ax^  term  should  be  solved 
for.  The  step  size  may  be  negative  or  positive  depending  on 
the  direction  of  travel. 

By  using  the  tangent  vector  corresponding  to  the  zero 
root,  successive  points  on  the  boundary  are  found.  At  each 
new  point  found,  the  error  in  the  approximation  is  computed 
using  Eq.  (29) .  By  using  the  criteria  concerning  step  size 
adjustment,  the  step  size  may  be  altered  as  additional  points 
on  the  boundary  are  determined. 

Since  a  stability  region  may  be  infinitely  large,  limits 
on  the  design  parameters  are  imposed  to  restrict  the  stability 
region  search  to  a  finite  area  in  the  x^,xj  plane.  If  either 
parameter  exceeds  a  limit  as  the  boundary  is  traversed,  then 
the  mapping  along  that  direction  is  discontinued.  A  new 
point  on  the  boundary  is  found  by  searching  along  the  limit 
line.  By  imposing  limits  on  x^  and  Xj,  the  resulting  stabil¬ 
ity  region  is  confined  to  lie  within  a  rectangle  formed  by 
the  limit  lines. 

Example  II 

To  illustrate  the  method  of  boundary  mapping,  consider 
the  following  example: 


>  A2  +  F^X)  X  +  Fq  (X)  =0 

(62) 

F^X)  =  X1  -  X2 

(63) 

Fq (X)  =  0.25 (x1  +  x2) 

(64) 

X  =  (x1#x2) 

(65) 
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The  base  point  is  given  as 

P°  =  (xj,x°)  =  (3,2)  (66) 

Let  x^  represent  the  axis  of  abscissas  and  x2  represent  the 
axis  of  ordinates.  Assume  that  the  first  point  on  the 
boundary  (P^)  has  been  found  using  the  root  locus  technique 
and  is  given  as 

P1  =  (xj,x°)  =  (2,2)  (67) 

At  p\  the  roots  of  G  are 

Xj  =  0  +  j  (68) 

X2  =  0  -  j  (69) 

The  approximating  equation,  Eq.  (21) ,  becomes 


AX1  =-OG/8x1^/8G/3X)  ^Ax1-(3G/3x2^/3G/8X)  Ax2  (70) 


ax2  ^OG/axj^^/aG/ax)  ^Ax^OG/ax^aG/ax)  ax2  (7i) 


After  evaluating  the  partial  derivatives  at  P1,  the  above 
equations  become 


AX 


1 


(0.25  +  i) 

- 


(0.25  -  j) 
2j 


(72) 


ax2  = 


IQ-2;, -  j)  Ax  +  (0.25  +  j) 

2j  flXl  +  Jj  ax2 


(73) 


Since  both  X^  and  X2  have  a  zero  real 
equation  may  be  used  to  determine  P  . 


part  at  P  ,  either 
By  using  the  AX^ 
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equation  and  rewriting,  the  following  equations  are 
obtained: 


Re  (AXX) 

=  -  0 . 5Ax^  + 

0 . 5 Ax 

(74) 

Im(AX^) 

=  0.125Ax1  + 

0.125Ax2 

(75) 

Requiring  that  Re ( A  X , )  =  0,  Eq.  (74)  represents  a  vector 
tangent  to  the  boundary  at  P  .  Since  the  slope  of  the 
tangent  vector  is  one,  let  AX2  equal  the  step  size  of  0.1 
and  solve  for  Ax^.  This  may  be  summarized  as  follows: 


Re(AX^)  =  -  0.5Ax^  +  0.5Ax.,  =  0 

(76) 

Ax^  =  Ax2  =  0.1 

(77) 

This  defines 

a  new  point  on  the  boundary  as 

P2  =  (x2,x2)  =  (2. 1,2.1) 

(78) 

2 

The  approximated  roots  at  P  are  determined  by 

t2  =  t1  +  at1 

(79) 

where  AX  is 

approximated 

determined  using  the  AX  given  by  Eq. 
■+2 

roots  at  P  are  then 

(77).  The 

\1  =  0  +  1. 025 j 

(80) 

X 2  -  0  -  1. 025  j 

(81) 

-*2 

The  error  estimate  at  P  may  be  computed  as  follows: 
error  =  X2  X^  -  Fq (X) 

error  =  (1. 025)  (1. 025)  -  1.05  =  0  (84) 
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Beginning  at  P  ,  additional  points  may  be  found  by  repeating 
this  process.  If  the  computed  error  estimate  at  any  point  is 
determined  to  be  larger  than  0.1,  then  G  is  re-solved  to 
provide  a  new  starting  point. 

Figure  3  is  the  stability  region  computed  for  this  exam¬ 
ple  using  the  method  of  tangent  vectors.  At  the  point  (5,5), 
the  search  along  the  upper  boundary  was  terminated  and  a  new 
point  on  the  boundary  was  found  at  (5,-5)  by  searching  along 
the  limit  line,  x^  =  5.  Table  3  summarizes  the  computa¬ 
tions  used  in  obtaining  Fig.  3. 


Figure  3.  Stability  region  map  for  Example  II. 
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Table  3.  Summary  of  Computations  used  in  Obtaining  Figure  3 


X1 

x2 

Tangent 
Vector ( X  ^ ) 

Tangent 

Vector (X2) 

Approxi¬ 

mate 

Roots 

Actual 

Roots 

2.0 

2.0 

(-0.5  4x^+0. 5 

Ax2) 

(-0.5  AXj+0.5  Ax2) 

- 

0  ±  1  j 

2.1 

2.1 

(-0.5  4x.+0.51 

Ax2) 

(-0. 514x^0. 51Ax2) 

0  -1. 025j 

0  ±  1 . 024 j 

2.2 

2.2 

(-O.SOIAXj+O.SOIAXj) 

(-0.5  Ax1+0 . 5  Ax2) 

0  ±1.05 

j 

0  ±  1.048 j 

2.4 

2.4 

(-0.5  Ax^+0. 5 

Ax2) 

(-0.5  AXj+0.5  Ax2) 

0  ±1.1 

j 

0  ±  1.09  j 

2.6 

2.6 

(-0.5  iXj+O. 5 

Ax2) 

(-0.5  4x^+0. 5  Ax2) 

0  ±1.15 

j 

0  ±  1.14  j 

3.0 

3.0 

(-0.5  AXj+0. 5 

Ax2) 

(-0.5  4x^+0 . 5  Ax2) 

0  ±1.24 

j 

0  ±  1.22  j 

3.4 

3.4 

(-0.5  AXj+0.5 

Ax2) 

(-0.5  Ax^+0.5  Ax2) 

0  ±1.32 

j 

0  ±  1.3  j 

4.2 

4.2 

(-0.5  AXj+0.5 

Ax2) 

(-0.5  AXj+0.5  Ax2) 

0  ±1.47 

j 

0  ±  1.48  j 

5.0 

5.0 

- 

- 

0  ±1.61 

j 

0  ±  1.58  j 

5.0 

-5.0 

(-AXj-AXj) 

(-0.984x^1. 03Ax2> 

- 

0,-10.0 

4.9 

-4.9 

(-AXj-AXj) 

(-0. 984x1+l. 034x2> 

0,-9. 8 

0,-  9.8 

4.8 

-4.8 

(-&Xj-Ax2) 

(-0.984x^1.014x2) 

0,-9. 6 

0,-  9.6 

4.6 

-4.6 

(-4x1-4x2) 

(-0.984x^1.024x2) 

0,-9. 2 

0,-  9.2 

4.4 

-4.4 

(-Ax^AXj) 

(-0.974x^1.034X2) 

0,-8. 8 

0,-  8.8 

4.0 

-4.0 

(-4x.-4x2) 

(-0.974x^+1.034x2) 

0 ,-8.0 

0,-  8.0 

3.6 

-3.6 

(-4x1-4x2) 

(-0.974x^1.034x2) 

0,-7. 2 

0,-  7.2 

2.6 

-2.6 

(-AXj-Ax2) 

(-0.964x^1.044X2) 

0,-5. 6 

0,-  5.6 

1.8 

-1.8 

(-Ax^-Ax2) 

(-0.934x^1.074X2) 

0,-3. 6 

0,-  3.6 

1.4 

-1.4 

(-Ax1-Ax2) 

(-0.914x^1.094X2) 

0,-2. 8 

0,-  2.8 

1.0 

-1.0 

(-Ax^-Ax2) 

(-0. 884x^+l. 124x2) 

0  ,-2 . 0 

0,-  2.0 

0.6 

-0.6 

(-AXj^-AXjj) 

(-0.794x^1.214X2) 

0,-1. 2 

0,-  1.2 

0.2 

-0.2 

(-AXj^-AXj) 

(-0.374xL+1.634x2) 

0,-0. 4 

0,-  0.2 

0.1 

-0.1 

(-Ax1-Ax2) 

(-0.254x^2.254x2) 

0,-0. 2 

0,-  0.2 

0.05 

0.05 

(-0.5  Ax1+0.5 

Ax2) 

(-0.5  AXj+0. 5  Ax2) 

- 

0  ±  0. 316 j 

0.15 

0.15 

(-0.5  4x^+0. 5 

Ax2) 

(-0.5  Ax^+0.5  4Xj) 

0  ±0 . 394 j 

0  ±  0.34  j 

0.25 

0.25 

(-0.5  4x^+0. 5 

Ax2) 

(-0.5  4Xj+0.5  Ax2) 

0  ±0.46 

j 

0  ±  0.41  j 

1.15 

1.15 

(-0.5  AXj+0.5 

Ax2) 

(-0.5  Ax^+0.5  4Xj) 

0  ±  0.75  j 
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4.3  REGION  CHARACTERISTICS 

The  stability  region's  boundary  was  defined  to  be  the 
set  of  points  within  the  plane  where  at  least  one  root  of 
the  characteristic  equation  is  zero  and  all  other  roots  are 
negative.  Within  the  stability  region,  the  roots  behave  in 
a  manner  predictable  by  the  tangent  vectors.  Lines  of  con¬ 
stant  damping,  constant  A,  and  constant  frequency  may  be 
constructed.  For  a  parametric  study  of  a  system  containing 
two  design  parameters,  a  two-dimensional  stability  region 
map  completely  describes  the  effects  of  parameter  variations. 

To  illustrate  the  information  which  can  be  obtained 
concerning  the  stability  region,  consider  the  previous  exam¬ 
ple.  As  indicated  in  Table  3,  the  tangent  vectors  coincide 
along  the  upper  boundary.  This  could  be  expected  since  both 
roots  have  a  zero  real  part  along  this  part  of  the  boundary. 

On  the  lower  boundary,  however,  the  tangent  vectors  are  not 
equal,  and  only  the  vector  corresponding  to  the  zero  root  is 
used  to  map  the  boundary.  By  using  the  vector  corresponding 
to  the  non-zero  root,  other  points  within  the  region  may  be 
obtained  which  yield  no  change  in  the  non-zero  root.  For 
example,  at  the  point  (1,-1),  the  non-zero  root  is 

X2  =  -  2.0  (83) 

The  corresponding  tangent  vector  is 

Re(AX2)  =  -  0.88Ax1  +  1.12Ax2  (84) 

Requiring  that  Re(AA2)  equal  zero  and  that  Re(AA1)  be  nega¬ 
tive  will  yield  a  point  within  the  region  where  A2  remains 
-2.0.  By  continuing  to  use  the  tangent  vector  for  \2»  addi¬ 
tional  points  may  be  obtained  where  A  2  is  a  constant  -2.0. 
Similarly,  constant  A  lines  for  different  A's  may  be  constructed. 
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Figure  4  is  a  detailed  stability  map  for  Example  II 
showing  the  construction  of  several  constant  X  lines  through¬ 
out  the  region.  It  is  interesting  to  note  that  each  constant 
X  line  is  tangent  to  a  parabola  which  corresponds  to  zero 
frequency.  The  parabola  consists  of  multiple  roots  and 
serves  to  separate  the  stability  region  into  two  areas.  The 
area  bounded  above  by  the  region's  upper  boundary  and  below 
by  the  multiple  root  parabola  contains  points  where  the  roots 
are  complex.  The  area  bounded  above  by  the  multiple  root 
parabola  and  below  by  the  region's  lower  boundary  contains 
points  where  the  roots  are  real. 

The  shaded  area  shown  in  Fig.  4  contains  those  points 
within  the  region  where  the  roots  are  complex.  Within  the 
complex  zone,  lines  of  constant  damping  may  be  constructed 
using  the  tangent  vectors.  For  example,  at  the  initial  base 


point  (3,2)  the  roots  are 

X1  =  -  0.5  +  j  (85) 

X2  -  -  0.5  -  j  (86) 

The  corresponding  tangent  vectors  are 

AXX  =  (0 . 5+0 . 125 j ) Ax1  +  (-0 . 5-0 . 375 j ) Ax2  (87) 

AX2  =  (0.5-0.125j) Axx  +  (-0.5+0. 375j)Ax2  (88) 


Since  the  real  parts  of  X1  and  X2  are  equal,  either  equation 
may  be  used.  Using  Eq.  (87)  and  rewriting  give 

Re(AX^)  &  0.5Ax^  =  0.5Ax2  (89) 

Requiring  Eq.  (89)  to  be  zero  will  yield  other  points  within 
the  region  where  the  real  part  of  X^  and  X2  is  -0.5.  Figure  4 
shows  several  lines  of  constant  damping  which  were  constructed 
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using  the  tangent  vectors.  The  lines  of  constant  frequency 
were  constructed  using  a  tangent  vector  corresponding  to  the 
imaginary  part  of  Eq.  (87).  For  example,  at  the  point  (3,2), 

ImtAX^  =  0.125Ax1  -  0.375Ax2  (90) 

Requiring  Eq.  (90)  to  be  zero  will  yield  other  points  where 
the  complex  part  of  X  remains  constant. 

4.4  THREE-DIMENSIONAL  REGIONS 

Two-dimensional  stability  regions  were  determined  by 
permitting  only  two  design  parameters  to  vary.  Three- 
dimensional  maps  can  be  obtained  by  generating  a  series  of 
two-dimensional  maps  for  different  values  of  a  third  param¬ 
eter.  The  map  obtained  in  this  manner  consists  of  a  series 
of  contour  curves  and  is  just  a  repetition  of  the  two- 
dimensional  process. 

5.0  EFFECTS  OF  PARAMETER  VARIATIONS  ON 
AIRCRAFT  LATERAL  STABILITY 


5.1  INTRODUCTION 

The  equations  of  motion  governing  an  aircraft  may  be 
linearized  such  that  they  form  two  independent  sets  of  equa¬ 
tions  pertaining  to  longitudinal  and  lateral  motion  (Ref.  4). 
The  linearized  equations  are  ordinary  linear  homogeneous 
differential  equations  whose  coefficients  are  functions  of 
various  aerodynamic  and  inertia  factors.  As  mentioned  previ¬ 
ously,  the  solution  to  such  a  set  of  equations  is  exponential 
and  an  investigation  of  stability  can  be  performed  by  ana¬ 
lyzing  the  eigenvalues  of  the  coefficient  matrix. 
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A  computer  program  using  the  root  movement  technique 
has  been  written  to  investigate  the  effects  of  variations 
in  the  aerodynamic  factors  on  the  lateral  stability  charac¬ 
teristics.  The  characteristic  equation  which  is  used  (Ref.  5) 
is  in  body  axis  and  is  of  the  form: 


where 


G  =  X4+F3(X)X3+F2(X)X2+F1(X)X+F0(X) 

X  =  (CVCVCVCVCjt_ 'CVCn g'Cn  'Cn 

pprpprppr 


(91) 

(92) 


The  exact  form  of  the  characteristic  equation  used  is  given 
in  Appendix  A. 


5.2  ROOT  LOCUS 


The  effects  of  parameter  variations  of  any  element  in 

*> 

X  are  accomplished  using  the  modified  root  locus  technique 
described  in  Section  3.0.  In  the  root  locus  routine,  only 
one  parameter  is  to  be  investigated,  and  hence,  there  is 
only  one  non- zero  element  in  the  vector  (AX) .  The  sparseness 
of  the  AX  vector  enables  one  to  investigate  only  the  appro¬ 
priate  column  of  the  matrix  (A).  For  example,  if  the  third 
element  of  3  is  to  be  investigated,  then  A^  contains  zeros 
in  all  positions  except  the  third.  This  requires  that  only 
the  third  column  of  A  be  evaluated  for  such  an  investigation 
and  eliminates  a  considerable  number  of  useless  computations. 


The  computer  output  from  the  root  locus  investigation 
is  the  value  of  the  roots  at  prescribed  increments  of  the 
parameter.  If  the  initial  base  point  is  stable,  then  the 
values  of  the  parameter  and  the  roots  are  printed  at  the 
first  point  of  instability.  In  addition  to  the  tabulated 
data,  the  program  has  a  printer  plot  routine  which  provides 
a  graphical  display  of  the  movement  of  the  roots  in  the 
complex  plane. 
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In  addition  to  tabulating  and  plotting  the  roots  at 
specified  increments  of  the  design  parameter,  the  forms  of 
the  zero  equation  and  its  roots  are  tabulated.  The  zero 
equations  for  each  of  the  nine  aerodynamic  terms  are  obtained 
by  applying  Eq.  (24)  to  Eq.  (91) .  These  equations  are  given 
in  Appendix  A. 

5.3  STABILITY  REGION  MAPPING 

Two-dimensional  stability  region  maps  can  be  obtained 
using  the  procedure  described  in  Section  4.0.  As  mentioned 
previously,  a  point  on  the  boundary  is  defined  to  be  a  point 
in  the  x^,xj  plane  where  at  least  one  of  the  roots  of  the 
characteristic  equation  has  a  zero  real  part.  Corresponding 
to  each  zero  root,  there  exists  a  tangent  vector  of  the  form 

Re(AXk)  =  a1Axi  +  a2Ax..  (93) 

Requiring  Eq.  (93)  to  equal  zero  will  yield  another  point  in 
the  x^,Xj  plane  where  Ak  remains  zero.  In  Eq.  (93),  Ax^  and. 
Axj  are  unknowns  and  a-L,a2  are  computed  using  Eqs.  (58)  and 
(59).  Since  Eq.  (93)  has  two  unknowns,  it  is  necessary  to 
specify  one  and  then  solve  for  the  other.  The  quantity 
specified  is  referred  to  as  the  step  size.  To  determine 
which  Ax  term  in  Eq.  (93)  is  to  be  specified  and  which  is 
to  be  solved  for,  it  is  necessary  to  examine  the  possibil¬ 
ities  which  may  occur  at  a  boundary  point.  Figure  5  illus¬ 
trates  the  possibilities  which  exist.  In  addition  to  knowing 
the  tangent  vector  slope,  it  is  necessary  to  know  the 
required  travel  direction.  The  direction  of  travel  is 
always  such  that  the  stable  region  is  to  the  right  of  the 
directed  line  segment  (PP^) .  In  Fig.  5,  the  direction  of 
travel  is  shown  to  be  positive  for  each  possibility,  where 
positive  travel  is  defined  to  be  the  direction  of  increasing 
ordinate  values.  For  negative  travel,  the  stable  zone  would 
be  opposite  that  pictured. 
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Figure  5.  Tangent  vector  possibilities. 
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The  solution  to  Eq.  (93)  then  provides  a  new  point  (P  ) 
on  the  boundary.  At  iP-,  the  error  in  the  roots  is  determined 
using  Eq.  (29).  If  the  error  is  excessive,  then  Eq.  (91)  is 
re-solved  and  the  step  size  is  halved. 

If  only  one  root  has  a  zero  real  part  at  the  boundary 
point,  then  its  corresponding  tangent  vector  is  used  to 
determine  a  neighboring  point  where  the  root  remains  zero. 
However,  if  more  than  one  root  is  zero,  then  the  various 
possibilities  which  arise  must  be  analyzed  and  the  proper 
tangent  vector  selected.  Note  that,  for  a  pair  of  complex 
roots  with  zero  real  part,  the  corresponding  real  tangent 
vectors  are  equal,  and  hence,  there  is  only  one  equivalent 
tangent  vector. 

If  at  a  boundary  point  there  are  two  different  tangent 
vectors  corresponding  to  zero  roots,  then  the  possibilities 
may  be  divided  into  two  classifications  according  to  slope 
signs.  If  the  two  tangent  vectors  have  opposite  slope 
signs,  then  the  correct  vector  and  travel  direction  may  be 
determined  by  examining  the  magnitude  of  the  slopes.  The 
various  possibilities  which  exist  are  illustrated  in  Fig.  6. 
Similarly,  if  the  two  tangent  vectors  have  the  same  slope 
signs,  the  correct  vector  and  travel  direction  may  be  deter¬ 
mined  by  examining  the  slope  magnitudes.  Figure  7  illustrates 
the  various  possibilities  which  exist  for  this  case.  It  can 
be  seen  that  for  the  dual  direction  possibilities,  the  problem 
of  determining  the  correct  tangent  vector  corresponds  to 
selecting  the  vector  which  is  to  the  right  of  the  last  move¬ 
ment  on  the  boundary.  In  both  Figs.  6  and  7,  V  represents 
the  last  tangent  vector  traveled.  At  the  point  (P1) ,  the  two 
possible  vectors  are  denoted  and  V2.  The  shaded  area 
represents  the  stable  part  of  the  region. 
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Figure  6.  Dual  direction  possibilities,  opposite  slope  signs. 
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Figure  7.  Dual  direction  possibilities,  same  slope  signs. 
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For  an  arbitrary  shaped  region,  points  on  the  boundary 

may  be  found  where  the  tangent  vector  undergoes  a  change  in 

slope  signs.  At  such  corner  points,  it  is  necessary  to 

explore  the  region  about  the  point  to  determine  the  correct 

travel  direction.  Figure  8  illustrates  the  various  corner 

point  possibilities.  In  Fig.  8,  V  represents  the  last  vector 
-►n 

traveled,  P  represents  the  corner  point,  represents  the 
tangent  vector  corresponding  to  the  zero  root  at  point  ?n, 
and  P*  represents  an  exploratory  point.  At  P*,  the  charac¬ 
teristic  equation  is  solved  to  determine  system  stability. 

By  using  the  information  acquired  at  the  exploratory  point, 
the  correct  travel  direction  along  V1  is  computed. 

The  stability  region  boundary  was  defined  to  be  the  set 
of  points  in  the  Xi,Xj  plane  where  at  least  one  of  the  roots 
of  the  characteristic  equation  has  a  zero  real  part.  In 
actual  practice,  the  defining  roots  are  required  only  to  be 
near  zero.  Hence,  the  set  of  points  computed  as  the  boundary 
lies  within  a  band  about  the  true  boundary.  Figure  9  illus¬ 
trates  an  arbitrary  region  where  the  true  boundary  is  indi¬ 
cated  by  the  solid  curve.  The  broken  curves  indicate  the 
boundary  zone. 

It  is  possible  in  traversing  the  boundary  to  lose  contact 
with  the  boundary  zone.  For  example,  at  point  P^  in  Fig.  9, 
the  tangent  vector  is  denoted  by  V.  By  using  the  tangent 
vector,  a  point  (P*)  is  determined  which  is  found  to  lie 
outside  the  boundary  zone.  By  re-setting  to  point  P^  and 
decreasing  the  step  size,  a  point  (P^+1)  can  be  found  which 
does  lie  within  the  boundary  zone. 

Since  stability  regions  may  extend  to  infinity,  realis¬ 
tic  limits  must  be  imposed  on  the  region  which  is  to  be 
mapped.  If  a  limit  is  encountered,  the  boundary  search  along 
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Stable  at  P*  Not  stable  at  P*  Stable  at  P*  Not  stable  at  P* 

travel  -  travel  +  travel  +  travel  - 


Figure  8.  Corner  point  tangent  vector  possibilities. 
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Figure  9.  Arbitrary  stability  region  map. 


the  last  direction  is  discontinued  and  a  one-dimensional 
search  is  made  along  the  limit  line  to  determine  a  new  point 
on  the  boundary.  Figure  9  illustrates  the  procedure  which 
is  used  when  a  limit  is  encountered.  The  arrowheads  indicate 
the  direction  of  travel. 

The  boundary  mapping  is  terminated  when  the  perimeter 
of  the  region  is  traversed.  This  is  indicated  when  the 
distance  from  the  point  (?n)  and  the  starting  point  (i*1)  is 
less  than  a  prescribed  tolerance.  If  the  point  (?n)  lies 
within  the  circle  encompassing  the  starting  point  t?1) ,  the 
mapping  process  is  terminated. 
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5.4  PROGRAM  OUTPUTS 

In  the  example  presented,  the  following  base  point  data 
were  used  and  are  representative  of  a  typical  twin- jet 


fighter  aircraft  (Ref. 

6)  . 

r  =  -  0.5730 

V 

0.2460 

1.1060 

C0  =  -  0.1570 

e 

ct  - 

p 

-  0.2170 

0.3050 

C  =  0.0808 

c 

nP 

0.0 

c»r- 

'  -  0.606 

nB 

The  additional  aircraft  and  : 

flight  condition 

data  used 

m  =  1215.0 

b  = 

38.41 

S  = 

538.34 

Ix  =  29,950 

Xz  = 

169,538 

Ixz 

5,241 

P/PQ  53  0.4485 

<*«,  = 

87.8 

V  = 

406.0 

a  =  15  deg 

0  = 

15  deg 

At  the  base  point,  the  roots  of  the  characteristic 
equation  are: 

X1  =  -  0.106  X2  =  -  0.496 

X 3  =  -  0.18  +  1.71J  X4  =  -  0.18  -  1 . 71j 

and  are  referred  to  in  the  following  discussion  as  the  poles. 
The  roots  correspond  to  the  spiral  mode,  the  rolling  mode,  and 
the  Dutch  roll  motion.  The  root  locus  plots  shown  in  Fig.  10 
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Figure  10.  Sample  root  locus  plot  (Cj^)  obtained  using 
the  lateral  stability  computer  program. 


illustrate  the  effects  that  variations  in  several  stability 
derivatives  have  on  the  lateral  stability  characteristics. 

The  zero  equations  corresponding  to  each  of  the  aerodynamic 
parameters  are  determined  using  Eq.  (24)  and  are  summarized 
in  Appendix  A. 

Figure  11  illustrates  a  stability  region  map  obtained 

using  the  computer  program  for  the  aerodynamic  parameters 

(Cn  and  C/  ) .  The  base  point  is  denoted  and  the  first 
p  p  -».i 

point  obtained  on  the  boundary  is  denoted  P  .  The  linear 

character  of  the  boundary  is  a  result  of  the  coefficient 

functions  of  the  characteristic  equation  being  linear  in 

Cn  and  .  The  origin  represents  a  multiple  root  point. 

p  $ 

Figure  12  illustrates  constant  A  lines  which  have  been 
constructed  throughout  the  stability  region.  The  X  lines 
correspond  to  the  spiral  and  rolling  mode  roots  (X^  and  X2) . 

It  can  be  seen  that  the  right  part  of  the  boundary  represents 
the  part  of  the  region  where  the  spiral  mode  becomes  unstable. 
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Figure  11.  Sample  stability  region  map  obtained  using 
the  lateral  stability  computer  program. 


Figure  12.  Detailed  stability  map  illustrating 
the  spiral  and  roll  roots. 
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The  line  of  X  =  -0.37  represents  a  line  of  multiple  roots. 
All  along  this  line,  the  spiral  and  rolling  mode  roots  are 
equal  and  become  complex  in  the  area  shaded. 


Figure  13  illustrates  constant  frequency  and  damping 
lines  for  the  Dutch  roll  roots.  It  can  be  seen  that  the 
left  part  of  the  boundary  represents  the  points  in  the 


Figure  13.  Detailed  stability  map  illustrating  the  Dutch  roll  roots. 

6.0  CONCLUSIONS 

The  technique  presented  herein  demonstrates  a  logical 
method  for  conducting  a  parametric  study  for  systems  which 
can  be  described  by  linear,  constant-coefficient  differential 
equations.  A  digital  computer  program  has  been  written 
which  applies  the  technique  to  determining  how  simultaneous 
variations  in  aerodynamic  parameters  affect  the  transient 
response  characteristics  of  an  aircraft. 
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As  previously  mentioned,  the  method  and  techniques 
described  can  be  applied  to  an  n^  degree  characteristic 
equation.  A  linearized  five  degree-of-freedom  model  is 
currently  being  developed  to  investigate  variations  in 
aerodynamic  parameters.  The  five  degree-of-freedom  model 
simplifies  to  the  three  degree-of-freedom  model  previously 
described . 


REFERENCES 

1.  Kaplan,  W.  Ordinary  Differential  Equations.  Addison- 

Wesley  Publishing  Company,  Reading,  Massachusetts, 
1967. 

2.  Kaplan,  W.  Advanced  Caluclus.  Addison-Wesley  Publishing 

Company,  Reading,  Massachusetts,  1959. 

3.  Rosenbach,  J.  B.,  Whitman,  E.  A.,  Meserve,  B.  E.,  and 

Whitman,  P.  M.  College  Algebra.  Fourth  Edition. 

Ginn  and  Company,  Boston,  Massachusetts,  1958. 

4.  Etkin,  B.  Dynamics  of  Flight.  John  Wiley  and  Sons,  Inc., 

New  York,  1959. 

5.  Stapleford,  R.  L. ,  Johnston,  D.  E.,  Teper,  G.  L.  and 

Weir,  D.  H.  "Development  of  Satisfactory  Lateral- 
Directional  Handling  Qualities  in  the  Landing 
Approach."  National  Aeronautics  and  Space  Administra¬ 
tion  CR-239,  Washington,  D.  C.,  July  ]965. 

6.  Chambers,  J.  R. ,  and  Anglin,  E.  L.  "Analysis  of  Lateral- 

Directional  Stability  Characteristics  of  a  Twin- Jet 
Fighter  Airplane  at  High  Angles  of  Attack."  National 
Aeronautics  and  Space  Administration  TN-D-5361, 

Langley  Research  Center,  Langley  Station,  Hampton, 
Virginia,  August  1969. 


46 


AEDC-TR-75-46 


APPENDIX  A 

LATERAL  STABILITY  CHARACTERISTIC  EQUATION 

G  =  X4+F3(X) A3+F2(X)X2+F1(X)X+F0 (X) 

where 

F3  (2)  — V-Ci  d.p+C2Kp+Nr+C3Lr) 

F2(X)  =  C1tYv(Kr+C3Lr+Lp+C2N  )+(N6+C3LB)  (cosci-Yr) 
-(Lj+C2Ns)  (sina+Yp)  +  (Nr+C3Lr)  (L  +C2N 
-G (Np+C3Lp) <Lr+C2Nr) ] 

Fl (X)  =  C1(YVC1[ (Np+C3Lp) (Lr+C2Nr)-(Nr+C3Lr) <Lp+C2Np) ] 

+  (Ng+CjLg)  [Cj^Yj.  (Lp+C2N_ )  -C^COSQ  (L  tCjN I  ) 

-Ssine-C^fL^C^)  -C^sinot (Lr+C2Nr) ] 

+ (Lg+C2Ng ) [-YrC1 (Np+C3Lp) +C1cosa (Np+C3Lp) 

-ScosQ+C^p  {Nr+C3Lr)+C1sina  CNr+C3Lr)  ]  } 

Fq(X)  =  C^{2cos0[Nr+C3Lr)  (Lg+C2Ng)-(NB-C3Lg)  (Lr+C2Nr)] 

-|sin9 [(Lg+C2Ng)  (Np+C2Lp)-(Lp+C3Np)  (Ng+C3Lg)]}  . 
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The  design  parameter  zero  equations  are: 

3G/3x1  =-X3+C1tx9+C3x6+x5+C2x8n2+c^t<x8+C3x5)  (Xg+CjXg) 

" (X9+C3X6*  ^X5+C2X8*  ^  X 

3G/3x2  =s""^1  tx4"^'^'2x7^  t"  (x^+C^x^ )  (Xg+C2xg) 

+ (x4+C2x7) (Xg+C^Xg) ] X 

3G/3x3  =-C^ Ix^+CgX^] X  +C^ [ (x7+C3x4) (Xg+C2Xg) 

- (x4+C2x7) (xg+C3x5) ] X 

3G/3x4  =  C1[C3(cosa-x3)+(-sina-x2) ]X2+C1[C3{C1x3(x5+C2x8) 

-C^cosa  (x5+C2x8)  -Ssine-Cj^  (x6+C2xg) x2 
-C^ (Xg+C2Xg ) sina} 

+{-C1x3 (Xg+C3x5)+C1cosa(x8+C3x5)-^cos0 
+C1x2 (xg+C3Xg)+C1sina(xg+C3Xg) }]X 

[cos 0  ( 1“C2C3 )  xg-sin0  {  (Xg+C2Xg ) —  C3  (x,-+C3Xg )  }  ] 
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3G/3x5  — C^+Cj.  [x^+CjXg  (1-C2C3)  ]  A2+C2  [C1C3x1  (x6+C2x9) 
(Xg+C3x4)+ (x7+C3x4) (x3-cosa)+ (x4+C2x7) 

(-x3+cosa) C3] A+C2 [-2sin9{C2 (x4+C2x7) - (x7+C3x4) } ] 

3G/3Xg  =-C1C3A3+C1  [C3{x1+C1  (x5+C2Xg)  (Xg+C-jXg)  ]  A2 
2 

+0^  [x^  (x8+C3x5)  +  (x?+C3x4)  (-x2-sina) 

+ (x4+C2x7) (C3x2+C3sina) ] A 
+C2 [^cos0x7 (C2C3~1) ] 

2 

3G/3x7  =  ^ [cosa-x3-C2sina-C2x2] A  +C1 [{x3C1 (x5+C2Xg) 

-C1cosa (x5+C2Xg) -2sin0-C1 (Xg+C2Xg) x2 

-C^  (Xg+C2xg) sina}+C2{-C1x3 (Xg+C3x5)+C1cosa (Xg+C3x,-) 

-2cos0+C1 (xg+C3Xg) x2+C1 (xg+C3Xg) sina}] A 

+C22[cos0 (C2C3-1) Xg-sin0 {C2 (xg+C2x5) - (x5+C3Xg) } ] 
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3G/3Xg  =-C1C2X3+  [C1C2x1+C1C2  (C^g+C^gXg)  -Cl  (CjXg+C^Xg )  ]  X2 
+ [x1{C1 (C1x6+C1C2xg) “C1C2 (C1x9+C1C3x4) } 

+ (C1x7+C1C3x4) (C1C2x3-C1C2cosa)+ (C1x4+C1C2x?) 

(-CiX3+Cic°sa)  ]  X+  [-^sin0  {C.^  (C1x4+C1C2x7) 
-C1C3(C1x7+C1C3x4) }] 

3G/3xg  =-C3X3+  (C1x5+C1C2xg)  -C1C2  (C^g+Cj^Xg)  ]  X2 

+  [c1c2x1  (c1xg+c1c3x5 )  -c1x1  (c^+c^Xg ) 

+ (C1x7+C1C3x4) (-C1C2x2-C1C2sina) + (C1x4+C1C2x7) 
*C1X2+Clsinc^  ^  *+  ty0086"^  (C1x4+C1C2x7) 

-CiC2 (c^x7+c^c3x4 )  }  ]  • 
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% 

clp 

Clr 

Cn 

Cn6 

CnP 

Cnr 

CY 

Cy3 

CyP 

Cy 

xr 

Cl 


C3 

Fi(X) 

fy 

g 


NOMENCLATURE 

Wing  span,  ft 

Rolling  moment  coefficient,  Mjj/q^Sb 

Static  derivative  of  rolling  moment  due  to  sideslip, 

ac^/36 

Dynamic  damping  derivative  in  roll,  3C£/(3pb/2V) 

Dynamic  cross-derivative  of  rolling  moment  due  to 
yawing,  3C£/(3rb/2V) 

Yawing  moment  coefficient,  Mz/qmSb 

Static  derivative  of  yawing  moment  due  to  sideslip, 
3Cn/36 

Dynamic  cross-derivative  of  yawing  moment  due  to 
rolling,  3Cn/ { 3pb/2V) 

Dynamic  damping  derivative  in  yaw,  3Cn/(3rb/2V) 

Side  force  coefficient,  Fy/q^S 

Static  derivative  of  side  force  due  to  sideslip, 

3  Cy/ 3  B 

Dynamic  derivative  of  side  force  due  to  rolling, 

3CY/ ( 3pb/2V) 

Dynamic  derivative  of  side  force  due  to  yawing, 
3Cy/(3rb/2V) 

1 

2 


ixz^x 

Ixz/Iz 

Characteristic  equation  coefficient  function 
Side  force,  lb 
Gravitational  acceleration 
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Ix  Moment  of  inertia  about  longitudinal  body  axis,  slug-ft2 

Ixz  Product  of  inertia,  slug-ft2 

Iz  Moment  of  inertia  about  normal  body  axis,  slug-ft2 

Lp  Rolling  moment  stability  parameter,  (pSb2V/4lx) C£ 

P 

Lr  Rolling  moment  stability  parameter,  ( pSb2V/4lx) C|r 

Lg  Rolling  moment  stability  parameter,  (pSbV/2lx) C|q 

Mx  Rolling  moment,  ft- lb 

Mz  Yawing  moment,  ft- lb 

m  Airplane  mass,  slug 

Np  Yawing  moment  stability  parameter,  (pSb2V/4Iz)Cnp 

Nr  Yawing  moment  stability  parameter, ( pSb2V/4lz) Cnr 

Ng  Yawing  moment  stability  parameter,  ( pSbV/2Iz) Cn^ 

The  ith  pole  is  represented  by  a  point  in  the  complex 
plane  which  is  given  by  the  value  of  the  ith  eigen¬ 
value  at  the  base  point 

Pm  Design  parameter  values  at  the  mth  point, 

=  <x“  x” . x“) 

Design  parameter  values  at  the  mth  point  and  the 
value  of  the  kth  eigenvalue  at  that  point 

p  Rolling  velocity,  rad/sec 

Free-stream  dynamic  pressure,  lb/ft2 

r  Yawing  velocity,  rad/sec 

S  Wing  area,  ft2 

V  Velocity,  ft/sec 

X  Design  parameter  column  vector,  $  =  (xi,x2 , . . . ,x^) 

Incremental  design  parameter  column  vector 
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Yr 

Y8 

Zi 

a 

8 

0 

t 

tk 

At 

P 

Po 


Side  force  stability  parameter,  (pSb/4m)CYp 

Side  force  stability  parameter,  (pSb/4m)Cyr 

Side  force  stability  parameter,  (pSV/2m)CYg 

The  ith  zero  point 

Angle  of  attack,  deg  or  rad 

Angle  of  sideslip,  deg  or  rad 

Angle  of  pitch,  deg  or  rad 

Eigenvalue  column  vector,  1  =  (X1,\2»***r^n) 
Eigenvalues  at  the  kth  point,  Xk  =  (Xk, Xk, • • • , Xk) 
Incremental  eigenvalue  column  vector 
Mass  density  of  air,  slug/ft3 
Mass  density  of  air,  sea  level,  slug/ft3 
Radian  frequency 
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